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Abstract: In this paper wc make a connection between the study of crowd dynamics and 
concepts from thermodynamics. The basic component of our continuous model is the continuity 
equation for the density of people, which we complete by prescribing the velocity field. This 
velocity field includes a nonlocal term modelling interactions between individuals. To provide 
support for our modelling assumptions we wish to prove an inequality that resembles the Second 
Law of Thermodynamics. To this end wc define an entropy-like functional and show that its 
time derivative equals a positive dissipation term minus a corrector term. The latter term should 
be small for the time derivative of the entropy to be positive. In case of isotropic interactions 
the corrector term is absent. For the anisotropic case we support the claim that the corrector 
term is small by performing simulations for the corresponding particle system. In fact, this term 
is sufficiently small for the entropy still to increase. Moreover, we can show that the entropy 
converges in time towards a limit value. 

Keywords: Crowd dynamics. Walking, View angles, Thermodynamics, Entropy, Steady states. 
First-order systems. 



1. INTRODUCTION 

Studying the behaviour of people in a crowd is nowadays 
no longer simply an activity of psychologists and social 
scientists. During the last few decades, it became clear 
that understanding and predicting the dynamics of people 
moving around is extremely important when designing 
buildings and infrastructure, and when being responsible 
for the safety of visitors at a large-scale event. 

Physicists and mathematicians were the ones who started 
providing answers to the questions that were posed. Their 
approach was (and is) similar to the way in which they deal 
with the physical, non-living world around us. Inevitably, 
this yields models in which people are nearly treated as 
if they were non-living material obeying physical laws. 
Illustrative for this way of thinking is the social force 
model, see e.g. Helbing and Molnar (1995). We remark that 
more or less in parallel and in the same spirit, the study of 
vehicular traffic and of biological aggregations (e.g. animal 
groups) developed. 

The difficulty in these models is to determine what the 
constitutive relations are. In physics, the intuition for these 
would be provided by experiments. We do not want to 
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claim that doing experiments with pedestrians is impossi- 
ble, but in any case it is difficult. One of the issues is the 
reproducibility. 

In this paper we explore a different way of justifying 
our modelling assumptions. We are interested in showing 
that a crowd obeys an inequality like the Second Law of 
Thermodynamics (also called Clausius-Duhem inequality) 
and thus maximizes 'entropy'. If we succeed, then this in- 
creases the consistency and trustworthiness of our model. 
For more background on thermodynamic concepts and 
their context, the reader is referred e.g. to Miiller and 
Ruggeri (1998). Maximizing entropy sounds nearly as if 
individuals try to win a game or maximize profit. See also 
Hoogendoorn and Bovy (2003) in this respect. 

We first describe our first-order continuum model (Section 
2), in particular addressing the constitutive relation we 
provide to close the model. This is the definition of the 
velocity field. Next, we introduce in Section 3 a concept 
of (generalized) entropy. If all interactions are isotropic 
with respect to the direction in which another individual 
is perceived, then one can show theoretically that the 
time derivative of the entropy is positive (in fact, non- 
negative). More work is required in the anisotropic case. 
Our numerical illustration in Section 4 suggests however 
that even in that case the entropy inequality holds. This is 



the main conclusion drawn here. We close the paper with 
an outlook on future work. 



2. MODEL EQUATIONS 



Consider the continuity equation 

dp 



dt 



+ V • [pv] = 



(1) 



on X 1R+ [d e N+ fixed). The model we use is very much 
in the spirit of Cristiani ct al. (2011), be it that their model 
is formulated in a more general way. Note that we have a 
first-order model, since we prescribe our velocity directly; 
see Coscia and Canavesio (2008) for an exposition of first- 
order models versus second-order models. We assume that 
the velocity field is the sum of two contributions 

v:=Vd + Vs, (2) 

a desired velocity v^, and a social velocity Vg- The desired 
velocity is taken to be a constant, while 



Vs{x) 



x-y 



VW{\x-y\)p{v)dv. (3) 



\x-y\ \vd\ 

The term Vs models the interactions between individuals. 
Note that we often omit the explicit time dependence of 
our variables. Here, W : M+ — >■ R is the potential governing 
the interactions. We use the word 'potential' here, since 
first-order models can be viewed as overdamped limits of 
second-order models. In the latter, the WW corresponds to 
a (gencrahzed) force. As such, W is a potential. We refer 
to Mogilner et al. (2003), p. 360, for a derivation in terms 
of a particle system. 
Here, 

x-y 



VW{\x-y\) := W'{\x-y\)-. 



(4) 



The fimction g : [—1,1] — > [0,1] incorporates anisotropy 
in the model. We have in mind the anisotropy that arises 
due to the fact that people have front and back sides. It 
depends on the direction in which one person perceives 
other people, how much influence they have on his motion. 
We restrict ourselves to linear functions g; that is, linear 
I , which is (—1 times) the cosine of the angle 



\x-y\ ivdi' 

under which point x sees point y. See also Gulikers et al. 
(2012) for our previous investigations on the effect of g on 
the dynamics of the underlying particle system. 

For the sake of being complete, we mention that the 
work by Coscia and Canavesio (2008) eventually does not 
focus on the nonlocal dependence on p, like in (3). In the 
cases they investigate in detail, they only allow for local 
dependence on p and/or Vp. 

The choice of this velocity field is an ansatz. A way of 
trying to justify this modelling assumption, is proving 
that the model is somehow consistent with ideas from 
thermodynamics. This will be the main motivation for all 
(mainly formal) steps in the sequel. 

As the continuity equation expresses the concept of con- 
servation of mass, we define the total mass 

M:= [ p{x) dx. (5) 

Moreover, we define the centre of mass 



Xq 



x) dx. 



(6) 



and the velocity of the centre of mass (or: barycentric 
velocity) 

Finally, we define the velocity with respect to the velocity 
of the centre of mass 



v{x) := v{x) — Vq. 



3. CLAUSIUS-DUHEM-LIKE INEQUALITY 



(8) 



For deriving an entropy-like functional (and corresponding 
Clausius-Duhem-like inequality) we first define and inves- 
tigate the following dissipation function D : [0, T] M+ 



D{t) := / p{x)\v{x)f dx. 



(9) 



Here, T > is some fixed final time. Note that, since Vd is 
constant, we have that 



^Ud^ I p(.z)dz . 
M .Ld M 



=Vd + 



M 



r7 / Vg{z)p{z)dz 
Vs{z)p{z)dz, (10) 



so 



V{x) -Vq= Vs{x) - 4j f '"a{z)p{z) dz. (11) 

A4 J'^d 

The second term on the right-hand side of (11) is indepen- 
dent of X. For the dissipation, we now have that 

D{t) = / p{v — vo)-vdx 

= / pVa-vdx — — - / Vspdz- / pvdx 

= pVs- vdx, (12) 

since Jjg.dPv dx = 0, which follows from the definition of i). 
For the ease of notation, let us define 

for all ^ G M'' \ {0}. If we take D and replace x by y (and 

vice versa), we obtain 

m = 



Pi^) l^J J{x - y)VW{\x - y\)p{y) dyj ■ v{x) dx 

/ P{y) i g{y-x)VW{\y-x\)p{x)dx) ■v{y)dy 

Jw^ \Jr<' / 

Pix) / giy-x)VW{\y-x\)-viy)p{y)dydx, (14) 



by changing the order of integration in the last step. We 
conclude from (4) that 

VW{\y-x\) = -VW{\x-y\). (15) 

Thus, 



D{t) = 

- P{x) I 9{y-x)VW{\x-y\)-v{y)p{y)dydx. 

(16) 

A combination of (12) and (16) yields 

2D{t) = I p{x) I VW{\x-y\yV{x,y)p{y)dydx, (17) 
where 

V{x, y) ■=v{x)g{x - y) - v{y)g{y - x) 
={v{x) - v{y)) 



^g{x-y) + ]^g{y-x) 



+ {v{x)+v{y)) 



\g{x -y)- \g{y - x) 



=:{v{x) - v{y))gs{x - y) + (w(x) + v{y))g^{x - y). 

(18) 

Here, gg and are the symmetric and antisymmetric parts 
of g, respectively. Thus 



2D{t)= / p{x) / VW{\x-y\)- 
Jr<i Jr'' 

■ {v{x) - v{y))gs{x - y)p{y) dy dx 
+ / p{x) I VW{\x-y\)- 



■ {v{x) + v{y))ga.{x - y)p{y) dy dx 
=:ss{t) + s,{t). (19) 

Inspired by Carrillo and Moll (2009), we define the follow- 
ing (entropy-like) functional 

S{t):=l [ p{x) [ W{\x-y\rg{x-y)p{y)dydx. (20) 

As the function g is linear, we can relate S to Sg- We first 
note that if g is linear, then g^ is constant. Let us call this 

constant a. 

Lemma 1. If the function g is linear, then 

dS_ _ 
~di ~ 



t; I P{x) I VW{\x-y\)-{v{x)-v{y))p{y)dydx. 

^ JR<^ JMA 



(21) 



Proof. 

dS 1 [ dp{x) 



dt 2 .Ld dt 



- y\)g{x - y)p{y) dydx 



+ \l P{x) I W{\x-y\y9{x-y)^dydx. 

(22) 

Interchange the order of integration in the first term, and 
replace a; by y (and vice versa) in the second term. Since 
W{\y — x\) = W{\x — y\), we obtain 



dS 1 

dt 2 jp^d 

1 



P{y) / W{\x-y\)g{x-y) 

Jr'' 



dp{x) 
dt 



+ ^ / P{y) I W{\x-y\)~g{y-x) 



dxdy 

dp{x 



dxdy 



p{y) I W{\x-ymx-y)^dxdy 



a / P{y) W^x-yDW ■ {p{x)v{x))dxdy, 

(23) 

where, in the last step, we used the continuity equation and 

= (X. Now, apply integration by parts in the x variable 
(where we assume vanishing boundary terms) to obtain 

dS 

— =a 
dt 



'■ / p{y) / ^W{\x-y\)-v{x)p{x)dxdy 



p{y) VWilx-yl) ■v{x)p{x)dxdy 



p{x) yW{\y-x\)-v{y)p{y)dydx 
= ^ / P{x) / VW{\x-y\)-v{x)p{y)dydx 



~% \ Pi^) '^W{\x-y\)-v{y)p{y)dydx. 

^ JR<^ JR<>- 

(24) 

In the last step we interchanged the order of integration 
in the first term, and used (15) in the second term. 
All this yields 

dS _ 

Itt ~ 

^ P{x) I WW{\x-y\)-{v{x)~v{y))p{y)dydx, 

(25) 

and the desired result then easily follows from the obser- 
vation that 

v{x)-v(y) = {v(x) + vo)-{v{y) + vo) = v{x)-v{y). (26) 

□ 

Remark: One can relax the linearity condition on g and 
obtain the same result. It suffices to have g'{ri) = g'{—r]) 
for all rj G [—1, 1]. Then V^s = and the corresponding 
term after integration by parts vanishes. For general (dif- 
ferentiable) g, this term remains, and thus an extra term 
appears in (21) (and eventually in the entropy inequality 
we are deriving). These technicalities are however beyond 
the scope and aim of this paper. 

The following equation summarizes our ideas so far in a 
condensed form: 



Dit)=^-^+s.it), 



dt 



(27) 



or 



dSjt) 
dt 



D{t) - s^{t), 



(28) 

which is a more desirable form, as it leads us to an entropy- 
like inequality. We first observe (by its definition (9)) that 
D{t) ^ for all t. Moreover, we identify the special case 
g = 1, in which we have fully isotropic interactions. To 
see this, observe in (3) that for = 1 the interaction term 
in Vg is VVF(|x — y\), which is the gradient of a radially 
symmetric function. Note that in this case gg = g and 



ffa = 0. 

obtain 



Thus, Sa{t) = for all t. As a consequence, we 
dS{t) 



dt 



= D{t) > 0, 



(29) 



which is a Clausius-Duhem-type inequality. This is a spe- 
cial case of what was treated by Carrillo and AfoU (2009), 
although there the reduction by subtracting vq is not done. 
If W is bounded, we have an upper bound on S (which is 
uniform in time). This implies that S will tend to some 
limiting value as t — )• oo. 

Our main question is now whether similar conclusions 
can be drawn if g is not constant. This is the case in 
which anisotropy is present in the interactions. In other 
words, some directions have more influence than others. 
In general, will no longer be 0. However, if this term is 
small (compared to D), still 

dS{t) 



dt 



^0 



(30) 



holds, which is the inequality we are looking for. In the 
sequel, we test this hypothesis numerically for a specific 
particle system. 

4. NUMERICAL ILLUSTRATION OF THE 
ANISOTROPIC CASE 

For a numerical illustration and investigation of the ideas 
described before, we simulate the particle system (of size 
N) corresponding to the model in Section 2. In particular, 
we take d= 2, and choose W to be the Morse potential 



S/Ir- 



(31) 



Furthermore, we take 



9{V) :=^(l + <7)-^(l-a)r?. 



(32) 



Here a G [0, 1] is a parameter that is used to tune the 
amount of anisotropy. If cr = 1, this corresponds to the 
isotropic case. Note that a relates to the aforementioned 
a via Q!=(l + cr)/2. We test the relatively simple case of 
N = 25, and take only cr = 0.5 for the sake of brevity. The 
latter choice ensures that we probe the anisotropic case. 

Initially the particles (individuals) are distributed ran- 
domly over the unit square [0, 1]^. We first show results 
for one simulation run, and afterwards for a sequence of 
1000 runs. New initial conditions are generated in each 
run. 

We want to see whether our numerics comply with dS/dt ^ 
0. To this end, we show in Fig. 1 the evolution in time 
of dS/dt. We deduced this quantity in three different 
ways from our simulations: Sg, the time derivative of S 
(where S was calculated for the particle system and the 
numerical time derivative was computed afterwards), and 
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1. Plots of Sg, dS/dt (numerical derivative of S), and 
D — Sa as functions of time. Theory predicts that they 
should be the same, as is supported by the graphs. It 
is important to note that the curves are positive and 
decay to zero (fluctuations around zero are 0(0.1)). 
Results for a single simulation run. 



W{s) ■.= Cae- 

see Mogilner et al. (2003), p. 363, or D'Orsogna et al. 
(2006) and the references cited therein. Note that D'Orsogna 
et al. (2006) use this potential in a second-order model. We 
demand that the parameters are positive and obey Ir < la 
and Cr/lr > Ca/la- This makes sure that the interactions 
arc repulsive in the short range, and attractive in the long 
range. (NB: These conditions include the Case 4 mentioned 
in MogUner et al. (2003): Cr > Ca and la > Iv) 
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Fig. 2. Plot of as a function of time. The value decays 
to zero. Results for a single simulation run. 

D — Sa. Theoretically these are identical, and the plots 
show that this is also the case numerically. We will thus 
not bother about this issue any more in the sequel. The 
main conclusion we draw from Fig. 1 is that dS/dt is indeed 
positive and moreover, that it decays to zero. 

Next, we examine the behaviour in time of the dissipation 
D. By deflnition, D is positive sec (9) - but it turns out 
that D actually tends to zero as t increases. This is shown 
in Fig. 2. Note that the behaviour of D is very similar to 
that of dS/dt. 

We observe in Figs. 1 and 2 that there is a sharp transition 
for f w 0.1 where the graphs become zero. If D — Sa and D 
are (nearly) zero, then ,Sa must also be (nearly) zero. The 
hypothesis that Sa is small, is thus valid after t — 0.1. For 
the interval [0, 0.1], we investigate the size of Sa separately. 
In Fig. 3 we show the ratio Sg^/D on the interval [0, 0.1]. We 
take this ratio since we are particularly interested in the 
size of Sa with respect to the size ofD. Indeed the values are 





Fig. 3. Plot of Si,/D as a function of time on the interval 
[0,0.1]. The values are small. Results for a single 
simulation run. 
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Fig. 4. Plot of 5 as a function of time. The value ap- 
proaches a limit value. Results for a single simulation 
run. 

small, i.e. of order 0(0.1). We do not continue the graph 
after t = 0.1 for a simple reason. Since both and D are 
very small then, we would divide two small numbers to 

get the ratio Sg_/D. The corresponding outcome does not 
provide any useful information. 

The interaction potential W we have taken is bounded. 
The entropy S therefore has a finite upper bound. As 
dS/dt is positive (cf. Fig. 1), we expect S to approach a 
limit value. This is indeed the case, as is shown in Fig. 4. 

We check now whether the conclusions we drew from one 
simulation run also hold if we do multiple runs (where in 
each run we impose different initial conditions). We do 
this to be sure that we have not just been 'lucky' so far. 
In Fig. 5 we show two curves corresponding to D — Sg,. 
At each time instance, we show both the minimum and 
maximum over all simulation runs. These graphs provide 
further evidence that D — (and hence dS/dt) is positive 
- which follows from the minimum - and decays to zero. 
Moreover, the term Sg, is thus small compared to D. After 
t ~ 0.25 the deviation from zero of the two curves is of 
order 0(0.1). 



Fig. 5. Plot of D — Sa as a function of time. At every time 
instance both the maximum and the minimum are 
taken over 1000 simulation runs. The graphs support 
the claim that D — Sa,\s positive and decays to zero. 
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Fig. 6. Plot of 5 as a function of time. The value ap- 
proaches a limit value. At every time instance the av- 
erage is taken over 1000 simulation runs. Wc checked 
that in all runs the same limit value was attained. 

Taking into account Fig. 5, which shows dS/dt > 0, we also 
expect S to increase and eventually converge to a limit. In 
Fig. 6 we show the average of S over the 1000 runs, and 
we see that the average indeed increases towards a limit. A 
plot of maxfe |S'i(t) — Sk(t)\ against time (where the index 
k runs over all simulation instances, and 5*1 is hence S 
obtained from the first simulation run), decays to zero. 
From this it follows that the limit value is the same in 
each simulation run. The graph is omitted here. 

5. CONCLUSIONS AND OUTLOOK 

We start this concluding section by answering the ques- 
tion: What does it actually mean for our system if the 
dissipation D goes to zero? By the definition of D in (9), 
D — automatically implies that v = 0. This means that 
all material points have the same velocity (namely the 
velocity of the centre of mass Vo{t)). For the continuum 
model, this means that the density profile p does not 
change shape any more; one can show that p is con- 



served along characteristics x{t) defined by the equation 

dx{t)/dt = vo{t). The configuration of the system is then 

just convected. One should realize however that this does 

not mean that the density p is uniform. 

By using the definition of 5 in (20), one can show that it 

automatically follows from D = that S is constant in 

time. 

We remark that the above explanation relates our entropy 
also to the theory of Lyapunov functionals. 

The numerics presented in this paper suggest that even 
in the case of anisotropic interactions, the dynamics still 
obey a Clausius-Duhem-type inequality. This supports the 
idea that it is worthwhile to do more effort to prove this 
analytically. We are aware of the fact that this fact might 
only be true under certain technical conditions, which were 
satisfied (more or less 'by accident') in our numerics. One 
of the aims of further theoretical investigations is to iden- 
tify these conditions, and make them as sharp as possible. 
Whereas the isotropic case is nearly a trivial exercise, it 
is clear from our (unsuccessful) attempts up to now that 
it will not be an easy task to prove the inequality in the 
anisotropic case. 

We stress here that proving an entropy inequality is not a 
goal in itself. It provides support for the 'thermodynamic 
consistency' of our model, and as such tells us that our 
ansatz for the velocity field is admissible, and hence, our 
modelling is not completely wrong. 

To get a better understanding of what exactly to prove, 
and how to proceed towards these proofs, we suggest to 
do first some more numerical experiments. In particular, 
we would pay attention to the following aspects: 

• One can do the same simulations, but change the 

value of a. Note that ,ga(-T — y) is proportional to 
(1 — (t), and this factor thus appears in Sa- The ampli- 
tude of this term most likely increases automatically 
with decreasing a {a 4- 0).^ Following this line of 
argument, it might be true that dS/dt is only positive 
for (7 within a certain distance from 1. This would lead 
to a condition on a; possibly a condition in which a 
is combined with other quantities. 

• Increase the number of particles N. One might guess 
that in a certain scaling, one can obtain the contin- 
uum model treated in this paper from the correspond- 
ing particle system that was used in the simulation 
section. The bigger N, the closer we then hope the 
result to be to the results of the continuum model. 
In this paper, a relatively small N was used, to have 
a system of ODEs that can still be handled easily. 
After all, the numerics in this paper only intend to 
illustrate our ideas and confirm our conjectures. 

• It would be interesting to see whether our results 
depend on the precise choice of interaction potential 
W. The simplest way would be to have different 
parameters {Cr,Ca,lrJa)- However, it is much more 
interesting to try another type of potential (possibly 
with a singularity around the origin). If the repulsive 

^ One should however be careful here in drawing this conclusion. 
Taking a closer to also changes the dynamics which, in turn, might 
cause a change in the integral term over p in Sa. This change could 
possibly counterbalance the change in 1 — cr in such a way that our 
claim on the size of Sa does not hold. 



behaviour around the origin is e.g. of the type ^ 1/r 
(cf. Coulomb interactions), we lose the trivial upper 
bound on S that followed from ||W^||oo- Numerics 
should then provide insight in whether the entropy 
still increases towards a limit value. Also, one could 
especially look at potentials that only model a zone 
of repulsion (that is, no attraction zone). This is 
interesting, since for repulsive interactions mass will 
'diffuse' and for each x fixed p{x,t) as t ^ oo. 
This implies that no steady states are to be expected. 
To what extent this destroys our entropy inequality 
is yet to be investigated. 
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